COMMENT ⊗   VALID 00007 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	BEGIN "EDGE_TEST"
C00003 00003
C00006 00004	IF EQU(ASK_STRING(" SKIP FALSE POSITIVES "),"Y") THEN GO TO AVG_DIFF
C00008 00005	for a range of conf values, find the corresponding 1% and 5%
C00011 00006	find the mean and std dev of s0sq/(noise_std↑2)
C00013 00007	now ask user which values he wants for diff and conf
C00016 ENDMK
C⊗;
BEGIN "EDGE_TEST"
	comment
		find best choice of conf and diff for Hueckel edge
		operator based on chosen false positive rate for
		no signal.  For a given false positive rate, find
		choice which requires the least height step to
		achieve a 50% false negative rate.  This uses
		steps generated by a Monte Carlo procedure;

DEFINE SAFE_ARRAY=1;	comment this makes all arrays safe;
REQUIRE "VISION.2[SYS,TOB]" SOURCE_FILE; comment vision language package;
REQUIRE EDGE_FILE;
REQUIRE EDGE_PRINT;
REQUIRE TV_VARIABLES;
REQUIRE "STEP4.SAI[FAK,TOB]" SOURCE_FILE;

comment
	conf=sqrt(s0sq/(s0sq+n0sq)) ie cosine
	diff=s0sq;

DEFINE NBHD=8;	comment the size of the array with fake intensity values;
LABEL FIRST_LOOP; ! entry for false positives and others;
LABEL FALSE_NEG; ! entry for false negatives;
LABEL AVG_DIFF;  ! entry for avg value of diff;
LABEL HISTO_LOOP; ! loop entry for histograms repeat with change scale;
DEFINE MAX_COUNT=2000; ! the max number of test samples generated;
INTEGER_ARRAY INTENSITY[1:NBHD*NBHD]; ! store fake intensities here;
REAL_ARRAY N0SQA[1:MAX_COUNT];	! store n0sq from edge oper here;
REAL_ARRAY S0SQA[1:MAX_COUNT];	! store s0sq from edge oper here;
INTEGER_ARRAY COUNT[0:100];	! histogram of scaled s0sq values;
SHORT INTEGER K,K1;
INTEGER	ENOUGH; ! how many samples (< max_count);
SHORT REAL STEP_MEAN,STEP_AMPLITUDE,NOISE_STD;
	comment the mean,amplitude,noise of fake steps;
SHORT REAL NSQ; ! scaling factor which depends on noise squared;
SHORT REAL NOISE_SCALE;	! scaling factor for histograms;
SHORT REAL THETA,RHO; ! angle,perp dist of edge from origin;
SHORT INTEGER BYTES;	! number of bytes per word;
SHORT REAL LAMBDA; ! wavelength (used to get phase near an edge);
SHORT INTEGER N1,N5,KN1,KN5; comment 1% and 5% level counts and indices;

ENOUGH←MAX_COUNT;
MODIFY(ENOUGH); ! ask user if he wants another value;
IF ENOUGH>MAX_COUNT THEN ENOUGH←MAX_COUNT;

comment now init the edge operator;
CONF←.05;DIFF←.05; ! choose low values to reject nothing;
BITS←6;
MODIFY(BITS);
LINLEN←WORDS(NBHD,BITS);
TVWORD←GIOWD(INTENSITY);
LSIDE←1;RSIDE←NBHD;
FLINE←1;LLINE←NBHD;
EJINIT;
DIFF←0;CONF←.05;

comment set up values for fake edges with zero height;
NOISE_STD←5.0; ! amplitude of gaussian noise;
MODIFY(NOISE_STD);
STEP_MEAN←(2↑BITS-1)*.5; ! offset to get positive intensities;
STEP_AMPLITUDE←0; ! step height for fake edges;
NOISE_SCALE←20.0;

IF EQU(ASK_STRING(" SKIP FALSE POSITIVES "),"Y") THEN GO TO AVG_DIFF;
FIRST_LOOP:
MODIFY(STEP_AMPLITUDE);

OUTSTR(CRLF&"FALSE POSITIVES ");
FOR K1←1 STEP 20 UNTIL ENOUGH DO BEGIN "ACCUM1"
FOR K←K1 STEP 1 UNTIL K1+20 DO BEGIN "ACCUM" ! report every 20 events;
	SHORT INTEGER N;
	LABEL CALLFK;
CALLFK:	FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
	STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
		comment generate an array with zero edge,random;
	EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
	N0SQA[K]←N0SQ;	! save the noise squared;
	S0SQA[K]←S0SQ;	! and save the signal squared;
	END "ACCUM";
PRINT(K);
END "ACCUM1";
comment for a range of conf values, find the corresponding 1% and 5%
false positives (above threshold with no edge signal) for fake
edges whose noise sq and signal sq have been saved;
N1←ENOUGH%100;
MODIFY(N1);
TCRLF;PRINT(N1);
HISTO_LOOP:
MODIFY(NOISE_SCALE);
NSQ←NOISE_SCALE/(NOISE_STD↑2);
FOR CONF←.95 STEP -.05 UNTIL .40 DO BEGIN "CONF_LOOP"	! loop over conf;
	SHORT REAL CONFSQ;
	SHORT INTEGER N;
	CONFSQ←CONF↑2;
	comment now form histograms over scaled s0sq values;
	FOR K←0 STEP 1 UNTIL 100 DO COUNT[K]←0; ! zero the histogram;
	FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ACCUM"
	SHORT REAL S;
	S←S0SQA[K];
	IF S/(S+N0SQA[K])>CONFSQ THEN BEGIN ! if above conf;
		N←S*NSQ+.5;
		IF N>100 THEN N←100;
		COUNT[N]←COUNT[N]+1;
		END;
	END "ACCUM";
OUTSTR(CRLF&" OUTPUT HISTOGRAM");
PRINT(CONF);
	INCHWL;
IFC 0 THENC
comment print histogram;
FOR K←0 STEP 1 UNTIL 100 DO BEGIN
	TCRLF;PRINT(K);OUTSTR(" "&CVS(COUNT[K]));
	END;
ENDC;

comment now find the 1% and 5% levels of the histogram tails;
N←0;
KN1←100;
FOR K←100 STEP -1 UNTIL 0 DO BEGIN
	N←N+COUNT[K];
	IF N<N1 THEN KN1←K;
	END;
OUTSTR(CRLF&"DIFF CUTOFF LIMITS "&CVG(N1/ENOUGH));
OUTSTR(CRLF&CVG(KN1/NSQ));
OUTSTR(CRLF&"DIFF CUTOFF LIMITS IN UNITS OF SQUARE STD DEV");
OUTSTR(CRLF&CVG(KN1/NOISE_SCALE));
PRINT(KN1);
END "CONF_LOOP";
IF EQU(ASK_STRING(CRLF&"REPEAT HISTOGRAMS "),"Y") THEN GO TO HISTO_LOOP;
IF EQU(ASK_STRING(CRLF&"REPEAT VALUES"),"Y") THEN GO TO FIRST_LOOP;
comment find the mean and std dev of s0sq/(noise_std↑2);

AVG_DIFF:
BEGIN "AVG_SIGNAL"
INTEGER K;
SHORT REAL SUM_S0SQ,SQ_S0SQ;
SHORT REAL STD_S0SQ; ! std dev of s0sq;
OUTSTR(CRLF&"MEAN VALUE OF DIFF AND STD DEV");
SUM_S0SQ←SQ_S0SQ←0;
TCRLF;
MODIFY(STEP_AMPLITUDE);
MODIFY(ENOUGH);
FOR K←1 STEP 20 UNTIL ENOUGH DO BEGIN "AVG_LOOP"
INTEGER K1;
PRINT(K);
FOR K1←1 STEP 1 UNTIL 20 DO BEGIN "A1_LOOP"
	FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
	STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
		comment generate an array with edge;
	EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
SUM_S0SQ←SUM_S0SQ+S0SQ;
SQ_S0SQ←SQ_S0SQ+S0SQ*S0SQ;
END "A1_LOOP";
END "AVG_LOOP";
STD_S0SQ←SQRT((SQ_S0SQ-SUM_S0SQ)/ENOUGH);
OUTSTR(CRLF&"AVERAGE VALUE OF DIFF/SQUARED AMPLITUDE"&CRLF);
OUTSTR(CVG(SUM_S0SQ/(ENOUGH*STEP_AMPLITUDE↑2)));
OUTSTR(CRLF&"AVERAGE VALUE OF STD DEV OF DIFF/SQUARED AMPLITUDE "&CRLF);
OUTSTR(CVG(STD_S0SQ/(STEP_AMPLITUDE↑2)));
NOISE_STD←0;
	FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
	STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
		comment generate an array with edge;
	EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
PRINT(NOISE_STD);PRINT(S0SQ);PRINT(STEP_AMPLITUDE);
END "AVG_SIGNAL";
comment now ask user which values he wants for diff and conf
start with low values of step height, low enough that fewer
than 50% of the fake edges are detected.  If too few, then
double step height.  If too many, we have bracketed step
height.  We can thus halve the step and use a binary search
to zero in on the right value;

FALSE_NEG:
BEGIN "EDGE_SEARCH" comment search for edge height with 50% detection;
REAL SEARCH_CONTROL;	comment either 2 or 1/2,
			start below, double value,bracket,binary search;
REAL STEP0,DSTEP;
LABEL CHANGE,EXIT,CHANGE1;
INTEGER K,J;	! loop count variables;


comment ask for initial values;
TCRLF;TCRLF;
STEP0←0;MODIFY(STEP0);
MODIFY(NOISE_STD);
NOISE_SCALE←1.0;
MODIFY(NOISE_SCALE);
DIFF←(NOISE_STD↑2)*NOISE_SCALE;
PRINT(DIFF);
CONF←0.90;MODIFY(CONF);

STEP_AMPLITUDE←STEP0;	COMMENT CHECK FALSE NEGATIVES;
SEARCH_CONTROL←2.0;
DSTEP←0.1*NOISE_STD;
THETA←PI/2;
RHO←0;
CHANGE1:
FOR J←1 STEP 1 UNTIL 10 DO BEGIN "TST1"
INTEGER TOTAL, POSITIVE; ! count of total tries and positives;
STEP_AMPLITUDE←STEP0+DSTEP;
TCRLF;MODIFY(STEP_AMPLITUDE);
MODIFY(DSTEP);
POSITIVE←0;
TOTAL←0;
MODIFY(ENOUGH);
FOR K←1 STEP 20 UNTIL ENOUGH DO BEGIN "TST2"
INTEGER K1;
PRINT(K);
FOR K1←1 STEP 1 UNTIL 20 DO BEGIN "TST3"
	FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
	STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISY_STEP,INTENSITY);
	IF EJLI(NBHD/2,NBHD/2,0,0)∧¬BCOMP THEN POSITIVE←POSITIVE+1;
	TOTAL←TOTAL+1;
	END "TST3";
	END "TST2";
PRINT(TOTAL);PRINT(POSITIVE);
	IF POSITIVE*2>TOTAL THEN SEARCH_CONTROL←0.5
		ELSE  STEP0←STEP0+DSTEP;
		DSTEP←SEARCH_CONTROL*DSTEP;
	END "TST1";

CHANGE:
	MODIFY(THETA);
	MODIFY(RHO);
	LAMBDA←1.0;
	GO TO CHANGE1;
	OUTSTR("ANOTHER VALUE-TYPE Y"&CRLF);	
	IF INCHWL≠"Y" THEN GO TO EXIT;
	MODIFY(RHO);
	MODIFY(THETA);
GO TO CHANGE;
EXIT:
END "EDGE_SEARCH";
END "EDGE_TEST";